//===================================
// G.Tatishvili
// PNNL
// mumupi0pi0_hbook-analysis code. Version 1.0.0
//===================================

#include <cstdio>
#include <cstdlib>
#include <iostream>

#include "CLHEP/Vector/LorentzVector.h"
#include "helix/Helix.h"

#include "event/BelleEvent.h"
#include "tuple/BelleTupleManager.h"

#include "basf/module.h"
#include "basf/module_descr.h"

#include "panther/panther.h"
#include BELLETDF_H
#include HEPEVT_H
#include MDST_H
#include MDST_OBS_H
#include EVTCLS_H
#include EVTVTX_H
#include LEVEL4_H
#include TRG_H

#include "ip/IpProfile.h"
#include "kid/atc_pid.h"
#include "eid/eid.h"
#include "mdst/mdst.h"
#include "mdst/Muid_mdst.h"
#include "mdst/findKs.h"

#include "kfitter/kvertexfitter.h"
#include "kfitter/kmassfitter.h"
#include "kfitter/kmassvertexfitter.h"
#include "kfitter/kmakemother.h"

#include "particle/Particle.h"
#include "particle/utility.h"
#include "particle/gammac.h"

#include "benergy/BeamEnergy.h"

#if defined(BELLE_NAMESPACE)
namespace Belle {
#endif

extern "C" {
  float prob_(const float&, const int&);
}


#define E_MASS   0.000511
#define MU_MASS  0.105658
#define PI_MASS  0.139570
#define K_MASS   0.493677
#define P_MASS   0.938272

#define LGNDR2(X) ((3.*X*X-1.)/2.)
#define LGNDR3(X) ((5.*X*X*X-3.*X)/2.)
#define LGNDR4(X) ((35.*X*X*X*X-30.*X*X+3.)/8.)
#define LGNDR5(X) ((63.*X*X*X*X*X-70.*X*X*X+15.*X)/8.)
#define LGNDR6(X) ((231.*X*X*X*X*X*X-315.*X*X*X*X+105.*X*X-5.)/16.)
#define LGNDR7(X) ((429.*X*X*X*X*X*X*X-693.*X*X*X*X*X+315.*X*X*X-35.*X)/16.)
#define LGNDR8(X) ((6435.*X*X*X*X*X*X*X*X-12012.*X*X*X*X*X*X+6930.*X*X*X*X-1260.*X*X+35)/128.)

  //fix below OR remove boost from HepParticles
  /*
  double E_HER = BeamEnergy::E_HER();
  double E_LER = BeamEnergy::E_LER();
  double cross_angle = BeamEnergy::Cross_angle();

  HepLorentzVector p4_her(sqrt(E_HER*E_HER-E_MASS*E_MASS)*sin(cross_angle),
                          0.0,
                          +sqrt(E_HER*E_HER-E_MASS*E_MASS)*cos(cross_angle),
                          E_HER);
  //p4_her.rotateY((float) cross_angle);
  HepLorentzVector p4_ler(0.0,
                          0.0,
                          -sqrt(E_LER*E_LER-E_MASS*E_MASS),
                          E_LER);
  HepLorentzVector P_BEAM = (p4_her + p4_ler);
  
  */

  HepPoint3D   IP(0,0,0);
  HepSymMatrix IPerr(3,0);

  atc_pid kid(3,1,5,3,2);
  atc_pid prid(3,1,5,4,2);

    void vertexFit(Particle &p, Particle &mu1, Particle &mu2)
    {
      kvertexfitter kvf;
      addTrack2fit(kvf,mu1); // add "K-" to kfitter.
      addTrack2fit(kvf,mu2); // add "PI+" to kfitter.
      if (IpProfile::usable()) addTube2fit(kvf);

      unsigned err = kvf.fit(); // do "fitting".
      if(err == 0){ // success.
	//    dynamic_cast<UserInfo&>(p.userInfo()).cl(kvf.cl());
	//dynamic_cast<UserInfo&>(p.userInfo()).chisq(kvf.chisq());
	//dynamic_cast<UserInfo&>(p.userInfo()).ndf(kvf.dgf());
	//	makeMother(kvf,p); // calculate "D0" using the fitting result.
	      std::cout << "kvf.vertex() "  << kvf.vertex() <<std::endl;
	      std::cout << "kvf.errVertex() "  << kvf.errVertex() <<std::endl;
	      std::cout << "kvf.chisq() "  << kvf.chisq() <<std::endl;

      }else{ // fail.
	//dynamic_cast<UserInfo&>(p.userInfo()).cl(-1.0);
	HepPoint3D vtx(999.,999.,999.);
	HepSymMatrix errVtx(3,0);
	p.momentum().decayVertex(vtx,errVtx);
      }
    }




  class HepParticle {
  public:

    HepParticle() {charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;mode=-1;charge=0;
      child[0]=child[1]=child[2]=child[3]=NULL;}

    HepParticle(Mdst_charged& trk, double mass, int type=-1) { // Constructer for charged particles
      if (type>=0 && type<=4) 
	{
	  Mdst_trk_fit& trkFit = trk.trk().mhyp(type);
	  
	  HepVector a(5,0);
	  a[0] = trkFit.helix(0);
	  a[1] = trkFit.helix(1);
	  a[2] = trkFit.helix(2);
	  a[3] = trkFit.helix(3);
	  a[4] = trkFit.helix(4);
	  
	  HepPoint3D pivot(trkFit.pivot(0), trkFit.pivot(1), trkFit.pivot(2));
	  Helix ltrk(pivot, a);
	  //	  ltrk.pivot(HepPoint3D(0.,0.,0.));
	  ltrk.pivot(IP);
	  //std::cout << "IP in particle declare" << IP << std::endl;
	  P4 = ltrk.momentum(0., mass);
	}
      else 
	{
	  HepLorentzVector P ( trk.px(), trk.py(), trk.pz(),
		      sqrt ( trk.px()*trk.px() + trk.py()*trk.py() +
			     trk.pz()*trk.pz() + mass*mass ) );
	  P4    = P;
	}

      gamma=NULL;pi0=NULL;vee2=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      charged = &trk;
      charge = (int)trk.charge();
      mode=-1;

      //            P4_boost = P4;
      //    P4_boost.boost(-P_BEAM.boostVector());
    }

    HepParticle(Mdst_gamma& gam) { // Constructer for gamma
      HepLorentzVector P ( gam.px(), gam.py(), gam.pz(),
		  sqrt ( gam.px()*gam.px() + gam.py()*gam.py() +
			 gam.pz()*gam.pz()) );
      P4    = P;
      charged=NULL;pi0=NULL;vee2=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      mode=-1;
      gamma = &gam;
      charge = 0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());

      //      std::cout << "P_BEAM 1 " << P_BEAM << std::endl;

    }

    HepParticle(Mdst_pi0& mdstpi0) { // Constructer for pi0
      HepLorentzVector P ( mdstpi0.px(), mdstpi0.py(), mdstpi0.pz(),
		  mdstpi0.energy() );
      P4  = P;
      charged=NULL;gamma=NULL;vee2=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      mode=0;
      pi0 = &mdstpi0;
      charge = 0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(Mdst_vee2& mdstvee2) { // Constructer for Ks
      HepLorentzVector P ( mdstvee2.px(), mdstvee2.py(), mdstvee2.pz(),
		  mdstvee2.energy() );
      P4  = P;
      charged=NULL;gamma=NULL;pi0=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      mode=0;
      vee2 = &mdstvee2;
      charge = 0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(HepParticle& ch1, HepParticle& ch2) { // 2-body Constructer
      P4=ch1.P4+ch2.P4;
      charge=ch1.charge+ch2.charge;
      child[0]=&ch1;
      child[1]=&ch2;
      child[2]=NULL;
      child[3]=NULL;
      charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;
      mode=0;


      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(HepParticle& ch1, HepParticle& ch2,
		HepParticle& ch3) {                   // 3-body Constructer
      P4=ch1.P4+ch2.P4+ch3.P4;
      charge=ch1.charge+ch2.charge+ch3.charge;
      child[0]=&ch1;
      child[1]=&ch2;
      child[2]=&ch3;
      child[3]=NULL;
      charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;
      mode=0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(HepParticle& ch1, HepParticle& ch2,
		HepParticle& ch3, HepParticle& ch4) { // 4-body Constructer
      P4=ch1.P4+ch2.P4+ch3.P4+ch4.P4;
      charge=ch1.charge+ch2.charge+ch3.charge+ch4.charge;
      child[0]=&ch1;
      child[1]=&ch2;
      child[2]=&ch3;
      child[3]=&ch4;
      charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;
      mode=0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }

  ~HepParticle() {};

    HepParticle & operator = (const HepParticle &p) { // Assigment
      P4 = p.P4;
      P4_boost = p.P4_boost;
      child[0] = p.child[0];
      child[1] = p.child[1];
      child[2] = p.child[2];
      child[3] = p.child[3];
      charged = p.charged;
      gamma = p.gamma;
      pi0 = p.pi0;
      vee2 = p.vee2;
      vertex = p.vertex;
      charge = p.charge;
      mode = p.mode;
    }
    

    void GetImpactParameters(double *dr, double *dz, int *rphi, int *z, int *ax, int *st, int t, HepPoint3D IPfrom = IP) {

      *dr = 1000.0;
      *dz = 1000.0;
      *rphi = 0;
      *z = 0;
      *ax = 0;
      *st = 0;

      if(charged->trk()){
	Mdst_trk_fit& trkFit = charged->trk().mhyp(t);
	//	std::cout << "SVD hitsrp " << trkFit.nhits(3) << std::endl;
	//std::cout << "SVD hitsz " << trkFit.nhits(4) << std::endl;
	if(trkFit){
	  HepVector a(5,0);
	  a[0] = trkFit.helix(0);
	  a[1] = trkFit.helix(1);
	  a[2] = trkFit.helix(2);
	  a[3] = trkFit.helix(3);
	  a[4] = trkFit.helix(4);

	  HepPoint3D pivot(trkFit.pivot(0), trkFit.pivot(1), trkFit.pivot(2));
	  Helix ltrk(pivot, a);

	  ltrk.pivot(IPfrom);

	  *dr = ltrk.dr();
	  *dz = ltrk.dz();
	  *rphi = trkFit.nhits(3);
	  *z = trkFit.nhits(4);
	  *ax = trkFit.nhits(0);
	  *st = trkFit.nhits(1);

	}
      }
    }
     
  // General informations
  HepLorentzVector P4;    		    	    	                // Four momentum in Lab-frame 
  HepLorentzVector P4_boost;		    	    	                // Four momentum in CM-frame
  HepParticle *child[4];                                        // Children
  Mdst_charged *charged;	    	    	                // Mdst banks for Pi,K
  Mdst_gamma *gamma;		    	    	                // Mdst banks for Gamma
  Mdst_pi0 *pi0;		    	    	                // Mdst banks for pi0
  Mdst_vee2 *vee2;  	    	    	    	                // Mdst banks for Vee2
  HepPoint3D vertex;                                            // vertex point
  int mode; 	    	    	    	    	                // mode code
  int charge;                                                   // charge of the particle
  };

  // Class definition
  class mumupi0pi0_hbook : public Module {
  public:
    mumupi0pi0_hbook ( void ){ eventcount=savecount=0;};
    ~mumupi0pi0_hbook ( void ){};
    void init ( int * ){}
    void term ( void ){
      printf("Final - %d/%d\n",savecount,eventcount); 
    };
    void disp_stat ( const char* ) {} ;
    void hist_def ( void );
    void event ( BelleEvent*, int* );
    void begin_run ( BelleEvent*, int * ){    
      eid::init_data();

      //get Beam Energy info
      BeamEnergy::begin_run();


      // get IP profile data from database
      IpProfile::begin_run();

      // dump IP profile data to STDOUT
      IpProfile::dump();    
    };
    void end_run ( BelleEvent*, int * ){};
    void other ( int*, BelleEvent*, int* ){}; 

private:
    BelleTuple *m_tpl, *m_tplpp;
    BelleTuple *m_new, *m_pi0, *m_2pi0, *m_pi0g, *m_2pi0g;
    BelleHistogram *h_hbmass,*h_hbprimemass;
    BelleHistogram *h_etabmass,*h_etabprimemass;
    BelleHistogram *h_Y5Smass,*h_Y4Smass,*h_Y3Smass,*h_Y2Smass,*h_Y1Smass;
    BelleHistogram *h_vtx, *h_vtxsv, *h_vtxprob;
    BelleHistogram *h_ntrk, *h_ntrksv, *h_dimu;

    void fwjet2(int NTrk, float *PTrk, float EBeam, 
		float *H0, float *R1, float *R2, float *R3, float *R4);


   
  int MCTruthInfo(HepParticle& Rec);
  int eventcount,savecount;
};

  //  extern float PROB(float chisq, int dof);

  extern "C" Module_descr *mdcl_mumupi0pi0_hbook ()
  {
    mumupi0pi0_hbook *module = new mumupi0pi0_hbook;
    Module_descr *dscr = new Module_descr ( "mumupi0pi0_hbook", module );
    IpProfile::define_global(dscr);
    return dscr;
  }
                                          
// Member Functions

// hist_def function
  void mumupi0pi0_hbook::hist_def ( void )
  {  
    extern BelleTupleManager *BASF_Histogram;
    BelleTupleManager& tm = *BASF_Histogram;
    BelleTupleManager& tmpp = *BASF_Histogram;
    BelleTupleManager& tmnew = *BASF_Histogram;
    BelleTupleManager& tmpi0 = *BASF_Histogram;
    BelleTupleManager& tm2pi0 = *BASF_Histogram;
    BelleTupleManager& tmpi0g = *BASF_Histogram;
    BelleTupleManager& tm2pi0g = *BASF_Histogram;
    
    /* Large nt for MC studies

    m_tpl = tm.ntuple ("Event", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2  \
pxpi1 pypi1 pzpi1 plabpi1 pxpi2 pypi2 pzpi2 plabpi2 pxg pyg \
pzg eg egboost recm2pi recmg2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2");

    m_tplpp = tmpp.ntuple ("Event pipi", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2 \
pxpi1 pypi1 pzpi1 plabpi1 pxpi2 pypi2 pzpi2 plabpi2 \
recm2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2");
*/

    //small ntupes for data
    m_tpl = tm.ntuple ("Event", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
vrp1 rz1 rp2 rz2 drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2 plabpi1 plabpi2 \
pzg eg egboost recm2pi recmg2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2 \
chisq chisqsv chisq2pi");

    m_tplpp = tmpp.ntuple ("Event pipi", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
rp1 rz1 rp2 rz2 drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2 plabpi1 plabpi2 \
recm2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2 \
chisq chisqsv chisq2pi");

    m_new = tmnew.ntuple ("Event new", "ntr npi0 ngam \
expno runno evtno energy prob2pi probev \
svd2pi angle2pi kprobpi1 kprobpi2 pprobpi1 pprobpi2 eid1 eid2 muid1 muid2 evtr2 recm2pi m2pi m2mu recmg2pi \
angleg eg egboost pqpi1 pqpi2 anglelabg dr1 dr2 dz1 dz2 rp1 rp2 rz1 rz2 ax1 ax2 st1 st2 pxpi1 pypi1 pzpi1 pxpi2 pypi2 pzpi2 recmpi1 recmpi2 npostr nnegtr mchrg npipi");

    m_pi0 = tmpi0.ntuple ("Event pi0", "ntr npi0 ngam \
expno runno evtno energy evtr2 recmpi npostr nnegtr mchrg pxpi1 pypi1 pzpi1");

    m_2pi0 = tm2pi0.ntuple ("Event 2pi0", "ntr npi0 ngam "
			    "expno runno evtno energy evtr2 recm2pi m2pi "
			    "recmpi1 recmpi2 npostr nnegtr mchrg npipi pxpi1 "
			    "pi1g1ecm pi1g2ecm pi2g1ecm pi2g2ecm "
			    "pypi1 pzpi1 pxpi2 pypi2 pzpi2 pi1g1e pi1g2e pi2g1e pi2g2e "
			    "pi1g1g2m p1g1p2g1 p1g1p2g2 p1g2p2g1 p1g2p2g2 pi2g1g2m "
			    "muid1 muid2 m2mu rcm2mup1 rcm2mup2 m2pi2mu recm2mu m2mupi1 m2mupi2");

//add eg egboost recmgpi0

    m_pi0g = tmpi0.ntuple ("Event pi0g", "ntr npi0 ngam \
expno runno evtno energy evtr2 recmpi npostr nnegtr mchrg eg egboost recmgpi0 pxpi1 pypi1 pzpi1");

//add eg egboost recmg2pi0

    m_2pi0g = tm2pi0g.ntuple ("Event 2pi0g", "ntr npi0 ngam \
expno runno evtno energy evtr2 recm2pi m2pi recmpi1 recmpi2 npostr nnegtr mchrg npipi eg egboost recmg2pi0 pxpi1 pypi1 pzpi1 pxpi2 pypi2 pzpi2");

    h_hbmass  = BASF_Histogram->histogram("hb MC mass", 100, 9.0, 10.0);
    h_etabmass  = BASF_Histogram->histogram("etab MC mass", 100,9.0, 10.0);
    h_hbprimemass  = BASF_Histogram->histogram("hb' MC mass", 100, 9.0, 10.0);
    h_etabprimemass  = BASF_Histogram->histogram("etab' MC mass", 100,9.0, 10.0);
    h_Y5Smass = BASF_Histogram->histogram("Y(5S) mass", 100, 9.0, 11.0);
    h_Y4Smass = BASF_Histogram->histogram("Y(4S) mass", 100, 9.0, 11.0);
    h_Y3Smass = BASF_Histogram->histogram("Y(3S) mass", 100, 9.0, 11.0);
    h_Y2Smass = BASF_Histogram->histogram("Y(2S) mass", 100, 9.0, 11.0);
    h_Y1Smass = BASF_Histogram->histogram("Y(1S) mass", 100, 9.0, 11.0);

    h_vtx = BASF_Histogram->histogram("vtx chisq", 102, -2, 100);
    h_vtxsv = BASF_Histogram->histogram("vtx chisq with sv", 102, -2, 100);
    h_ntrk = BASF_Histogram->histogram("vtx trk", 21, -0.5, 20.5);
    h_ntrksv = BASF_Histogram->histogram("vtx trk with sv", 21, -0.5, 20.5);
    h_vtxprob = BASF_Histogram->histogram("vtx prob with sv", 100, 0., 1.);

    h_dimu = BASF_Histogram->histogram("DiMuon Mass", 100, 9., 11.);

  }

// event function
void mumupi0pi0_hbook::event ( BelleEvent* evptr, int* status )
{
  // Get managers for MDST banks
  Mdst_charged_Manager& ChgMgr = Mdst_charged_Manager::get_manager();
  Mdst_gamma_Manager& GamMgr = Mdst_gamma_Manager::get_manager();
  Mdst_pi0_Manager& Pi0Mgr = Mdst_pi0_Manager::get_manager();
  Mdst_vee2_Manager& Vee2Mgr = Mdst_vee2_Manager::get_manager();
  Belle_event_Manager& EvtMgr = Belle_event_Manager::get_manager();
  Gen_hepevt_Manager& HepevtMgr = Gen_hepevt_Manager::get_manager();  
  Mdst_ecl_Manager& EclMgr = Mdst_ecl_Manager::get_manager();
  Mdst_ecl_aux_Manager& EclAuxMgr = Mdst_ecl_aux_Manager::get_manager();
  Evtcls_hadronic_flag_Manager& HadFlagMgr = Evtcls_hadronic_flag_Manager::get_manager();
  Evtcls_flag_Manager& FlagMgr = Evtcls_flag_Manager::get_manager();
  Evtcls_hadron_info_Manager& HadInfoMgr = Evtcls_hadron_info_Manager::get_manager();
  L4_summary_Manager& L4SumMgr = L4_summary_Manager::get_manager();
  Evtvtx_primary_vertex_Manager& VtxMgr = Evtvtx_primary_vertex_Manager::get_manager();
  Rectrg_summary3_Manager& TrgSumMgr = Rectrg_summary3_Manager::get_manager();

  *status=0;


  // set IP and error
  int IPUsable = 0;
  if(IpProfile::usable()){
    IP    = IpProfile::position(1);
    IPerr = IpProfile::position_err_b_life_smeared(1);
    IPUsable = 1;
  }else {
    IP    = HepPoint3D(0,0,0);
    IPerr = HepSymMatrix(3,0);
  }

  //--------------------------------

  /*  float e_cms = 10.8667;

  double E_HER_0 = 7.998213;
  double E_LER_0 = 3.499218;
  double E_BEAM  = e_cms/2;
  double E_HER = (E_BEAM/5.29)*E_HER_0;
  double E_LER = (E_BEAM/5.29)*E_LER_0;
  double cross_angle = 0.022;
  double MASS_E = 0.000511;
  */

  double E_BEAM = BeamEnergy::E_beam_corr();
  double E_HER = BeamEnergy::E_HER();
  double E_LER = BeamEnergy::E_LER();
  double cross_angle = BeamEnergy::Cross_angle();

  HepLorentzVector p4_her(0.0,
                          0.0,
                          +sqrt(E_HER*E_HER-E_MASS*E_MASS),
                          E_HER);
  p4_her.rotateY(cross_angle);
  HepLorentzVector p4_ler(0.0,
                          0.0,
                          -sqrt(E_LER*E_LER-E_MASS*E_MASS),
                          E_LER);
  HepLorentzVector P_BEAM = (p4_her + p4_ler);
  
  //  std::cout << "P_BEAM 2 " << P_BEAM << std::endl;

  //--------------------------------


  int Necl_barrel = 0;
  for(std::vector<Mdst_ecl>::iterator it = EclMgr.begin();
      it != EclMgr.end(); it++){
    Mdst_ecl& cls = *it;

    if(cls.quality() == 0 && cls.energy() > 0.1 && //E_min_cut 
       //if(cls.quality() == 0 && cls.energy() > 0.04 && //E_min_cut 
       cos(cls.theta()) > -0.7 && // ECLcos_min
       cos(cls.theta()) <  0.9 )  // ECLcos_max
      Necl_barrel++;
  }


  int evt_Ntrk = -1, evt_Ncls = -1;
  double evt_Psum, evt_Esum, evt_Evis, evt_Pz, evt_HeavyJetMass, evt_Thrust, evt_R2;
  evt_Psum = evt_Esum = evt_Evis = evt_Pz = evt_HeavyJetMass = evt_Thrust = evt_R2 = -1.;
  
  if (HadInfoMgr.count()>0) {
    evt_Ntrk = HadInfoMgr[0].Ntrk(); 
    evt_Ncls = HadInfoMgr[0].Ncls();
    evt_Psum = HadInfoMgr[0].Psum();
    evt_Esum = HadInfoMgr[0].Esum();
    evt_Evis = HadInfoMgr[0].Evis();
    evt_Pz   = HadInfoMgr[0].Pz();
    evt_HeavyJetMass = HadInfoMgr[0].HeavyJetMass();
    evt_Thrust    = HadInfoMgr[0].Thrust();
    evt_R2    = HadInfoMgr[0].R2();
  }

// Define Track list
  std::vector<HepParticle> TrkList,GamList,Pi0List, MupList, MumList;
  std::vector<Particle> PartTrkList,PartTrkpList,PartTrkmList;

  int Npos = 0;
  int Nneg = 0;
  
  HepLorentzVector P_CHG;

// Loop over charged tracks

//  std::vector<HepParticle>::iterator tit = TrkList.begin();
  for(std::vector<Mdst_charged>::iterator it = ChgMgr.begin(); it != ChgMgr.end(); ++it)
    {
      Mdst_charged& Charged = *it;


      //      std::cout<< "IP " << IP << std::endl;


      Particle Partpl(Charged,"PI+",IP); 
      Particle Partmi(Charged,"PI-",IP); 


      if( !good_charged(*it) ) continue;

      if (Charged.charge()>0)
	{ Npos++;            
	}
      if (Charged.charge()<0) 
	{Nneg++;     
	}

      //      HepParticle part(Charged,PI_MASS,2); 
      HepParticle part(Charged,MU_MASS,1); 

      Muid_mdst muID(Charged);
	

     
      

      //      Mdst_trk_fit &trk_fit = Charged.trk().mhyp(2);
      Mdst_trk_fit &trk_fit = Charged.trk().mhyp(1);

      //            if (Charged.trk().quality() == 0
	  //	  	  && trk_fit.nhits(3)>0
	  //&& trk_fit.nhits(4)>1
      //      	  )
	{
	  //	  std::cout << "PI_MASS="<<PI_MASS<<std::endl;
	  TrkList.push_back(part);
	  //	  if (Charged.charge()>0) {PartTrkList.push_back(Partpl);}
	  //if (Charged.charge()<0) {PartTrkList.push_back(Partmi);}
	  if (Charged.charge()>0) {PartTrkpList.push_back(Partpl);}
	  if (Charged.charge()<0) {PartTrkmList.push_back(Partmi);}


	  P_CHG = P_CHG + part.P4;


	  if (muID.Chi_2()>0){
	    double mu_id = muID.Muon_likelihood();
	      {
		if (Charged.charge()>0) {MupList.push_back(part);}
		if (Charged.charge()<0) {MumList.push_back(part);}
	      }
	  }
	}
    }



// Loop over gamma 
  for(std::vector<Mdst_gamma>::iterator it = GamMgr.begin(); it != GamMgr.end(); ++it)
    {    
      Mdst_gamma& Gamma = *it;           
      HepParticle part(Gamma);
      double cos_theta = part.P4.pz()/part.P4.rho();
      Mdst_ecl_aux &aux = EclAuxMgr( Panther_ID(Gamma.ecl().get_ID()) );
      double e9oe25 = aux.e9oe25();
	
      //      std::cout << "e9oe25=" << e9oe25 << std::endl;

      if ((Gamma.ecl().match() == 0 || Gamma.ecl().match() == 2)
	  //&&		 sqrt ( Gamma.px()*Gamma.px() +Gamma.py()*Gamma.py() +
	  //	   Gamma.pz()*Gamma.pz())> 0.5)
	  )
      {    
	if (e9oe25 < 0.95) continue;
	if (part.P4.e() < 0.04) continue;
// 	if (cos_theta > 0.8501 && part.P4.e() < 0.100) continue;
// 	if (-0.6643 < cos_theta && cos_theta < 0.8501 && part.P4.e() < 0.050) continue;
// 	if (-0.6643 < cos_theta && cos_theta < 0.8501 && part.P4.e() < 0.1) continue;
// 	if (cos_theta < -0.6643 && part.P4.e() < 0.100) continue;
	
	// v2 below
	/*	if (e9oe25 < 0.95) continue;
	if (cos_theta > 0.8501 && part.P4.e() < 0.080) continue;
	if (-0.6643 < cos_theta && cos_theta < 0.8501 && part.P4.e() < 0.040) continue;
	if (cos_theta < -0.6643 && part.P4.e() < 0.080) continue;
	*/

	//    if (cos_theta > 0.8501) continue; //not forward barrel


	HepParticle part(Gamma);
	GamList.push_back(part);
	}
    }


  int ntrk=TrkList.size()+GamList.size(),trkcount=0;
  float ptrk[ntrk*3],H0,R1,R2,R3,R4;

  for(std::vector<HepParticle>::iterator it = TrkList.begin();
      it != TrkList.end(); ++it){

    HepParticle& Charged = *it;

    ptrk[trkcount*3]=Charged.P4_boost.px();
    ptrk[trkcount*3+1]=Charged.P4_boost.py();
    ptrk[trkcount*3+2]=Charged.P4_boost.pz();
    trkcount++;
  }

  for(std::vector<HepParticle>::iterator it = GamList.begin();
      it != GamList.end(); ++it){

    HepParticle& Gamma = *it;

    ptrk[trkcount*3]=Gamma.P4_boost.px();
    ptrk[trkcount*3+1]=Gamma.P4_boost.py();
    ptrk[trkcount*3+2]=Gamma.P4_boost.pz();
    trkcount++;
  }

  fwjet2(ntrk,ptrk,E_BEAM,&H0,&R1,&R2,&R3,&R4);

  int pi0_count = 0;

  for(std::vector<Mdst_pi0>::iterator it = Pi0Mgr.begin();
      it != Pi0Mgr.end(); ++it){        
  
    Mdst_pi0& Pi0 = *it;    

    if (fabs(0.134-Pi0.mass())>=0.05*3) continue;

    HepParticle part(Pi0);  // Create a Pi0 as HepParticle  
    part.child[0] = part.child[1] = NULL;
    for(std::vector<HepParticle>::iterator it = GamList.begin();
	it != GamList.end(); ++it) {

      HepParticle& ga = *it;
      if (Pi0.gamma(0) == (*ga.gamma)) part.child[0] = &ga;
      if (Pi0.gamma(1) == (*ga.gamma)) part.child[1] = &ga;

    }
    
    
    if (part.child[0]!=NULL && part.child[1]!=NULL)
      {
	if (part.child[0]->P4.e() < part.child[1]->P4.e()) {
	  HepParticle* tmp = part.child[0];
	  part.child[0] = part.child[1];
	  part.child[1] = tmp;
	}
	  
// 	Mdst_gamma& gamma0=Pi0.gamma(0);
// 	Mdst_gamma& gamma1=Pi0.gamma(1);

	Mdst_gamma gamma0 = *(part.child[0]->gamma);
	Mdst_gamma gamma1 = *(part.child[1]->gamma);

	HepParticle gam0(gamma0);
	HepParticle gam1(gamma1);

	double cos_theta0 = gam0.P4.pz()/gam0.P4.rho();
	double cos_theta1 = gam1.P4.pz()/gam1.P4.rho();

	Mdst_ecl_aux &aux0 = EclAuxMgr( Panther_ID(gamma0.ecl().get_ID()) );
	double e9oe250 = aux0.e9oe25();

	Mdst_ecl_aux &aux1 = EclAuxMgr( Panther_ID(gamma1.ecl().get_ID()) );
	double e9oe251 = aux1.e9oe25();

	//	std::cout << "e9oe250,1 costheta0,1 " << e9oe250 << " " << e9oe251 << " " << cos_theta0 << " " << cos_theta1 << std::endl; 
	

	//	std::cout << "CRAIGB: gam0e = " << gam0.P4.e() << ", gam1e = " << gam1.P4.e() << std::endl;

	if (gam1.P4.e()<0.200){std::cout << "pi0 photon energy " << gam0.P4.e() << " " << gam1.P4.e() << std::endl;}
	if (e9oe250 < 0.95 && e9oe251 < 0.95) continue;
	//	if (cos_theta0 > 0.8501 && cos_theta1>0.851) continue;
	//	if (gam0.P4.e() < 0.1 &&  gam1.P4.e() < 0.1) continue;


	if (e9oe250 < 0.95) continue;
	if (gam0.P4.e() < 0.040) continue;
// 	if (cos_theta0 > 0.8501 && gam0.P4.e() < 0.100) continue;
// 	if (-0.6643 < cos_theta0 && cos_theta0 < 0.8501 && gam0.P4.e() < 0.050) continue;
// 	if (-0.6643 < cos_theta0 && cos_theta0 < 0.8501 && gam0.P4.e() < 0.1) continue;
// 	if (cos_theta0 < -0.6643 && gam0.P4.e() < 0.100) continue;

	if (e9oe251 < 0.95) continue;
	if (gam1.P4.e() < 0.040) continue;
// 	if (cos_theta1 > 0.8501 && gam1.P4.e() < 0.100) continue;
// 	if (-0.6643 < cos_theta1 && cos_theta1 < 0.8501 && gam1.P4.e() < 0.050) continue;
// 	if (-0.6643 < cos_theta1 && cos_theta1 < 0.8501 && gam1.P4.e() < 0.1) continue;
// 	if (cos_theta1 < -0.6643 && gam1.P4.e() < 0.100) continue;



	//for v2 
	/*	if (e9oe251 < 0.95) continue;
	if (gam1.P4.e() < 0.040) continue;
	if (cos_theta1 > 0.8501 && gam1.P4.e() < 0.080) continue;
	if (-0.6643 < cos_theta1 && cos_theta1 < 0.8501 && gam1.P4.e() < 0.040) continue;
	if (cos_theta1 < -0.6643 && gam1.P4.e() < 0.080) continue;

	if (e9oe250 < 0.95) continue;
	if (gam0.P4.e() < 0.040) continue;
	if (cos_theta0 > 0.8501 && gam0.P4.e() < 0.080) continue;
	if (-0.6643 < cos_theta0 && cos_theta0 < 0.8501 && gam0.P4.e() < 0.040) continue;
	if (cos_theta0 < -0.6643 && gam0.P4.e() < 0.080) continue;
	*/



	Pi0List.push_back(part);  
	pi0_count++;

	//	std::cout << "pi0count,energy " <<pi0_count<< " " << Pi0.energy() << std::endl;
      }
    // else 
    //   std::cout << "Never made it into Pi0 loop: " << part.child[0] << " " << part.child[1] << std::endl;
  }

  //  std::cout << "CRAIGB: " << Pi0Mgr.size() << " " << Pi0List.size() << std::endl;

  //make gamma list that is not included among pi0 photons

  std::vector<HepParticle> myGamList;




  for(std::vector<HepParticle>::iterator it0 = GamList.begin();
      it0 != GamList.end(); ++it0){
    
    HepParticle& Gamma0 = *it0;
    int mypi0_count = 0;

    
    for(std::vector<HepParticle>::iterator it1 = GamList.begin();
	it1 != GamList.end(); ++it1){
      
      HepParticle& Gamma1 = *it1;

      HepLorentzVector P_Pi0 = Gamma0.P4+Gamma1.P4;
      
      double Mass_Pi0 = P_Pi0.m();
      
      if (fabs(0.134-Mass_Pi0)<=0.05*3){mypi0_count++;}

      //      std::cout << "soft pi0 " << P_Pi0.rho() << std::endl;
      
    }

     if (mypi0_count==0)
      {
	HepParticle part(Gamma0);
	myGamList.push_back(part);
	/*if (Gamma0.P4.e()>=0.100&&Gamma0.P4.e()<0.150){
	  if (my_ngam[0]<9){ my_ngam[0]=my_ngam[0]+1;}}
	if (Gamma0.P4.e()>=0.150&&Gamma0.P4.e()<0.200){
	  if (my_ngam[1]<9){ my_ngam[1]=my_ngam[1]+1;}}
	if (Gamma0.P4.e()>=0.200&&Gamma0.P4.e()<0.250){
	  if (my_ngam[2]<9){ my_ngam[2]=my_ngam[2]+1;}}
	if (Gamma0.P4.e()>=0.250&&Gamma0.P4.e()<0.300){
	  if (my_ngam[3]<9){ my_ngam[3]=my_ngam[3]+1;}}
	if (Gamma0.P4.e()>=0.300&&Gamma0.P4.e()<0.350){
	  if (my_ngam[4]<9){ my_ngam[4]=my_ngam[4]+1;}}
	if (Gamma0.P4.e()>=0.350&&Gamma0.P4.e()<0.400){
	  if (my_ngam[5]<9){ my_ngam[5]=my_ngam[5]+1;}}
	if (Gamma0.P4.e()>=0.400&&Gamma0.P4.e()<0.450){
	  if (my_ngam[6]<9){ my_ngam[6]=my_ngam[6]+1;}}
	if (Gamma0.P4.e()>=0.450&&Gamma0.P4.e()<0.500){
	  if (my_ngam[7]<9){ my_ngam[7]=my_ngam[7]+1;}}
	if (Gamma0.P4.e()>=0.500&&Gamma0.P4.e()<0.550){
	  if (my_ngam[8]<9){ my_ngam[8]=my_ngam[8]+1;}}
	if (Gamma0.P4.e()>=0.550){
	  if (my_ngam[9]<9){ my_ngam[9]=my_ngam[9]+1;}}*/

	  	 
      }
  }
  //end attempt to make new gamma list

  
    

  Gen_hepevt_Manager &hepmgr = Gen_hepevt_Manager::get_manager();
  for (std::vector<Gen_hepevt>::iterator it = hepmgr.begin();
       it != hepmgr.end(); ++it) {
    
    if (abs(it->idhep()) == 10553){
      h_hbmass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 551){
      h_etabmass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 110553){
      h_hbprimemass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 100551){
      h_etabprimemass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 553){
      h_Y1Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 100553){
      h_Y2Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 200553){
      h_Y3Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 300553){
      h_Y4Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 9000553){
      h_Y5Smass->accumulate ( (float) it->M(), 1.0);
    }
  }

  //  if (Npos <1 || Nneg <1 ) return;

  if (Pi0List.size() == 0 ) return;

  // if ((MupList.size() > 0) && (MumList.size() > 0) && (Pi0List.size() > 1))
  //   std::cout << "CRAIGB: found an event! " << std::endl;

  //now for main loop


  	      int my_ngam[10];
	      int gam_bin = 0;
	      
	      for (gam_bin=0;gam_bin<10;++gam_bin){
		my_ngam[gam_bin]=0;
	      }


	      int my_npipi[10];
	      int pipi_bin = 0;
	      
	      for (pipi_bin=0;pipi_bin<10;++pipi_bin){
		my_npipi[pipi_bin]=0;
	      }



	      // find mu pair - then add one track, then add two tracks.

	      //vertex the two mu tracks
	      //	      kvertexfitter kv;

	      std::vector<Particle>::iterator tit3 = PartTrkpList.begin();		
	      for(std::vector<HepParticle>::iterator it3 = MupList.begin(); it3 != MupList.end(); ++it3)
		{
		  std::vector<Particle>::iterator tit4 = PartTrkmList.begin();
		  for(std::vector<HepParticle>::iterator it4 = MumList.begin(); it4 != MumList.end(); ++it4)
		    {
		      
		      kvertexfitter kv;

		      HepParticle& Mu1  = *it3;
		      HepParticle& Mu2  = *it4;

		      Muid_mdst Mu1_muid(*(Mu1.charged));
		      Muid_mdst Mu2_muid(*(Mu2.charged));

		      float mu_id1 = Mu1_muid.Muon_likelihood();
		      float mu_id2 = Mu2_muid.Muon_likelihood();
		      
		      HepParticle P2DIMU(Mu1,Mu2);

		      HepLorentzVector P_DIMU = Mu1.P4+Mu2.P4;
		      

		      Particle& tMu1  = *tit3;
		      Particle& tMu2  = *tit4;

		      addTrack2fit(kv, tMu1);
		      addTrack2fit(kv, tMu2);
		
		      if (IpProfile::usable()) addTube2fit(kv);

		      unsigned fitstatus = kv.fit(); // do "fitting".


  /*    	      std::cout << "trk # " << i << std::endl;
	      std::cout << "kvsv.fit() "  << kvsv.fit() <<std::endl;
	      std::cout << "kv.vertex() "  << kv.vertex() <<std::endl;
	      std::cout << "kv.errVertex() "  << kv.errVertex() <<std::endl;
	      std::cout << "kv.chisq() "  << kv.chisq() <<std::endl;

	      std::cout << "trk sv# " << isv << std::endl;
	      std::cout << "kvsv.fit() "  << kvsv.fit() <<std::endl;
	      std::cout << "kvsv.vertex() "  << kvsv.vertex() <<std::endl;
	      std::cout << "kvsv.errVertex() "  << kvsv.errVertex() <<std::endl;
	      std::cout << "kvsv.chisq() "  << kvsv.chisq() <<std::endl;
  */  

		      if(fitstatus == 0){ // success.
			Ptype ptype_ups("PSI");
			Particle p(P_DIMU, ptype_ups); 

		      if (P2DIMU.P4.m()>9)
			{


			vertexFit(p, tMu1, tMu2);

			//			makeMother(kv,p); // calculate "D0" using the fitting result.
			//			HepPoint3D vtx(kv.vertex());
			//HepSymMatrix errVtx(kv.errVertex());
			//			p.momentum().decayVertex(vtx,errVtx);
			//			std::cout << "p.momentum " << p.momentum().decayVertex(vtx,errVtx) << std::endl;


			//			Particle Mu1fit(tMu1, 1, kv.vertex());
			//Particle Mu2fit(tMu2, 1, kv.vertex());

			//			h_dimu->accumulate ( (float) p.mass(), 1.0);

			h_dimu->accumulate ( (float) P2DIMU.P4.m(), 1.0);

			}
		      }else{ // fail.
			//			dynamic_cast<UserInfo&>(p.userInfo()).cl(-1.0);
			HepPoint3D vtx(999.,999.,999.);
			HepSymMatrix errVtx(3,0);
			//			p.momentum().decayVertex(vtx,errVtx);


		      }


		      //		      std::cout << " particle p " << p.momentum() << std::endl;

		      //	
		      double dr1,dz1;
		      int nrp1, nz1;
		      int nax1, nst1;


		      /*		      std::cout << "before mu mass" << Mu1.P4.x() << std::endl;
		      std::cout << "mu mass" << Mu2.P4.x() << std::endl;

		      Mu1.GetImpactParameters(&dr1,&dz1,&nrp1,&nz1,&nax1,&nst1,0, IP);      
		      Mu2.GetImpactParameters(&dr1,&dz1,&nrp1,&nz1,&nax1,&nst1,0, IP);      

		      

		      std::cout << "after mu mass" << Mu1.P4.x() << std::endl;
		      std::cout << "mu mass" << Mu2.P4.x() << std::endl;
		      */

		      //		      h_dimu->accumulate ( (float) P2DIMU.P4.m(), 1.0);
		      //		      h_dimu->accumulate ( kv.momentum(), 1.0);

		      
		      //		      continue;
		      if (P2DIMU.P4.m()>9)
			{




  //loop over pi0

  for(std::vector<HepParticle>::iterator it1 = Pi0List.begin(); it1 != Pi0List.end(); ++it1)
    {
      HepParticle& Pi01  = *it1;

      HepLorentzVector P_RECp1  = P_BEAM-Pi01.P4;
      double m_rec1  = P_RECp1.m();


	  if (EvtMgr.count()>0) 
	    {
	      m_pi0->column("expno",EvtMgr[0].ExpNo());
	      m_pi0->column("runno",EvtMgr[0].RunNo());
	      m_pi0->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
	    }
	  m_pi0->column("ntr",     TrkList.size());
	  m_pi0->column("npi0",    Pi0List.size());
//	  m_pi0->column("ngam",    my_ngam_count);
	  m_pi0->column("ngam",    myGamList.size());
	  m_pi0->column("energy",E_BEAM);
	  m_pi0->column("evtr2",evt_R2);

	  m_pi0->column("pxpi1",   Pi01.P4.px());
	  m_pi0->column("pypi1",   Pi01.P4.py());
	  m_pi0->column("pzpi1",   Pi01.P4.pz());
	  
	  m_pi0->column("recmpi", m_rec1);


	  m_pi0->column("npostr",  Npos);
	  m_pi0->column("nnegtr",  Nneg);
	  
	  m_pi0->column("mchrg", (P_CHG).m());
	  
	  /*	Mdst_gamma& gamma0=Pi01.gamma(0);
	Mdst_gamma& gamma1=Pi01.gamma(1);

	HepParticle gam0(gamma0);
	HepParticle gam1(gamma1);

	double cos_theta0 = gam0.P4.pz()/gam0.P4.rho();
	double cos_theta1 = gam1.P4.pz()/gam1.P4.rho();

	Mdst_ecl_aux &aux0 = EclAuxMgr( Panther_ID(gamma0.ecl().get_ID()) );
	double e9oe250 = aux0.e9oe25();

	Mdst_ecl_aux &aux1 = EclAuxMgr( Panther_ID(gamma1.ecl().get_ID()) );
	double e9oe251 = aux1.e9oe25();


	  m_pi0->column("cosg0", cos_theta0 );
	  m_pi0->column("cosg1", cos_theta1 );
	  m_pi0->column("e9e250", e9oe250 );
	  m_pi0->column("e9e251",e9oe251);
	  m_pi0->column("eg0", gam0.P4.e());
	  m_pi0->column("eg1", gam1.P4.e());
	  */

	      if (m_rec1 > 9.2 && m_rec1 < 11.0){
	  m_pi0->dumpData();}
	  
	      for(std::vector<HepParticle>::iterator it2 = it1+1; it2 != Pi0List.end(); 
		  //for(std::vector<HepParticle>::iterator it2 = Pi0List.begin(); it2 != Pi0List.end(); 
	  ++it2)
	{

	  //	  if (it1==it2) continue;

	  HepParticle& Pi02  = *it2;

	  // Ensure Pi01 is the fast pi0
	  if (Pi01.P4.e() < Pi02.P4.e()) {
	    HepParticle tmp = Pi01;
	    Pi01 = Pi02;
	    Pi02 = tmp;
	  }

	  P_RECp1  = P_BEAM-Pi01.P4;
	  m_rec1  = P_RECp1.m();

	  HepParticle P2PI(Pi01,Pi02);
	  HepParticle P2MU2PI(P2DIMU,P2PI);

	  HepLorentzVector P_RECp2  = P_BEAM-Pi02.P4;
	  HepLorentzVector P_REC  = P_BEAM-Pi01.P4-Pi02.P4;
	  HepLorentzVector P_RECmu  = P_BEAM-Mu1.P4-Mu2.P4;
	  
	  double m_rec  = P_REC.m();
	  double m_recmu  = P_RECmu.m();
	  double m_rec2  = P_RECp2.m();

	  HepLorentzVector P_REC_com = P_REC;

	      if (P_REC.m()>=9.80&&P_REC.m()<9.85){
		if (my_npipi[0]<9){ my_npipi[0]=my_npipi[0]+1;}}
	      if (P_REC.m()>=9.85&&P_REC.m()<9.90){
		if (my_npipi[1]<9){ my_npipi[1]=my_npipi[1]+1;}}
	      if (P_REC.m()>=9.90&&P_REC.m()<9.95){
		    if (my_npipi[2]<9){ my_npipi[2]=my_npipi[2]+1;}}
	      if (P_REC.m()>=9.95&&P_REC.m()<10.0){
		if (my_npipi[3]<9){ my_npipi[3]=my_npipi[3]+1;}}
	      if (P_REC.m()>=10.0&&P_REC.m()<10.05){
		if (my_npipi[4]<9){ my_npipi[4]=my_npipi[4]+1;}}
	      if (P_REC.m()>=10.05&&P_REC.m()<10.15){
		if (my_npipi[5]<9){ my_npipi[5]=my_npipi[5]+1;}}
	      if (P_REC.m()>=10.15&&P_REC.m()<10.20){
		    if (my_npipi[6]<9){ my_npipi[6]=my_npipi[6]+1;}}
	      if (P_REC.m()>=10.20&&P_REC.m()<10.30){
		if (my_npipi[7]<9){ my_npipi[7]=my_npipi[7]+1;}}
	      if (P_REC.m()>=10.30&&P_REC.m()<10.35){
		if (my_npipi[8]<9){ my_npipi[8]=my_npipi[8]+1;}}
	      if (P_REC.m()>=10.35&&P_REC.m()<10.40){
		if (my_npipi[9]<9){ my_npipi[9]=my_npipi[9]+1;}}
	      

	      P_REC_com.boost(-P_BEAM.boostVector());

	      double hb_angle_com = P_REC_com.pz()/P_REC_com.rho();
	      double hb_angle = P_REC.pz()/P_REC.rho();

	      if (m_rec < 9.2 || m_rec > 11.0) continue;

	      if( evt_R2< 1.0 //<0.4
		  )
		{ //if cuts

	      savecount++;
                      
	      *status = 1;
		

	  int my_npipi_count=0;

	  for (pipi_bin=0;pipi_bin<10;++pipi_bin){
	    my_npipi_count=my_npipi_count+(int)pow(10,pipi_bin)*my_npipi[pipi_bin];
	  }
	      if (EvtMgr.count()>0) 
		{
		  m_2pi0->column("expno",EvtMgr[0].ExpNo());
		  m_2pi0->column("runno",EvtMgr[0].RunNo());
		  m_2pi0->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
		}
	      m_2pi0->column("ntr",     TrkList.size());
	      m_2pi0->column("npi0",    Pi0List.size());

//	      m_2pi0->column("ngam",    my_ngam_count);
	      m_2pi0->column("ngam",    myGamList.size());
	      m_2pi0->column("energy",E_BEAM);
	      m_2pi0->column("evtr2",evt_R2);

	      m_2pi0->column("pxpi1",   Pi01.P4.px());
	      m_2pi0->column("pypi1",   Pi01.P4.py());
	      m_2pi0->column("pzpi1",   Pi01.P4.pz());

	      m_2pi0->column("pxpi2",   Pi02.P4.px());
	      m_2pi0->column("pypi2",   Pi02.P4.py());
	      m_2pi0->column("pzpi2",   Pi02.P4.pz());

              m_2pi0->column("recm2pi", m_rec);
              m_2pi0->column("recm2mu", m_recmu);
              m_2pi0->column("m2pi",    P2PI.P4.m());

              m_2pi0->column("recmpi1", m_rec1);
              m_2pi0->column("recmpi2", m_rec2);

	      m_2pi0->column("npostr",  Npos);
	      m_2pi0->column("nnegtr",  Nneg);

	      m_2pi0->column("mchrg", (P_CHG).m());

	      m_2pi0->column("npipi",    my_npipi_count);


	      double mbc2;
	      mbc2 = (2*E_BEAM)*(2*E_BEAM)
		- (P2MU2PI.P4.px()-P_BEAM.px())*(P2MU2PI.P4.px()-P_BEAM.px())
		- (P2MU2PI.P4.py()-P_BEAM.py())*(P2MU2PI.P4.py()-P_BEAM.py())
		- (P2MU2PI.P4.pz()-P_BEAM.pz())*(P2MU2PI.P4.pz()-P_BEAM.pz());
	      
	      m_2pi0->column("muid1",mu_id1);
	      m_2pi0->column("muid2",mu_id2);

	      m_2pi0->column("pi1g1e", Pi01.child[0]->P4.e());
	      m_2pi0->column("pi1g2e", Pi01.child[1]->P4.e());
	      m_2pi0->column("pi2g1e", Pi02.child[0]->P4.e());
	      m_2pi0->column("pi2g2e", Pi02.child[1]->P4.e());

	      HepLorentzVector Pi01g1_cm(Pi01.child[0]->P4.px(), Pi01.child[0]->P4.py(),  
					 Pi01.child[0]->P4.pz(), Pi01.child[0]->P4.e());
	      HepLorentzVector Pi01g2_cm(Pi01.child[1]->P4.px(), Pi01.child[1]->P4.py(),  
					 Pi01.child[1]->P4.pz(), Pi01.child[1]->P4.e());
	      HepLorentzVector Pi02g1_cm(Pi02.child[0]->P4.px(), Pi02.child[0]->P4.py(),  
					 Pi02.child[0]->P4.pz(), Pi02.child[0]->P4.e());
	      HepLorentzVector Pi02g2_cm(Pi02.child[1]->P4.px(), Pi02.child[1]->P4.py(),  
					 Pi02.child[1]->P4.pz(), Pi02.child[1]->P4.e());

	      Pi01g1_cm.boost(-P_BEAM.boostVector());
	      Pi01g2_cm.boost(-P_BEAM.boostVector());
	      Pi02g1_cm.boost(-P_BEAM.boostVector());
	      Pi02g2_cm.boost(-P_BEAM.boostVector());
	       
	      m_2pi0->column("pi1g1ecm", Pi01g1_cm.e());
	      m_2pi0->column("pi1g2ecm", Pi01g2_cm.e());
	      m_2pi0->column("pi2g1ecm", Pi02g1_cm.e());
	      m_2pi0->column("pi2g2ecm", Pi02g2_cm.e());

	      m_2pi0->column("pi1g1g2m", (Pi01.child[0]->P4 + Pi01.child[1]->P4).m());
	      m_2pi0->column("p1g1p2g1", (Pi01.child[0]->P4 + Pi02.child[0]->P4).m());
	      m_2pi0->column("p1g1p2g2", (Pi01.child[0]->P4 + Pi02.child[1]->P4).m());
	      m_2pi0->column("p1g2p2g1", (Pi01.child[1]->P4 + Pi02.child[0]->P4).m());
	      m_2pi0->column("p1g2p2g2", (Pi01.child[1]->P4 + Pi02.child[1]->P4).m());
	      m_2pi0->column("pi2g1g2m", (Pi02.child[0]->P4 + Pi02.child[1]->P4).m());

	      HepLorentzVector P_REC3 = P_BEAM-Pi01.P4-Mu1.P4-Mu2.P4;
	      double m_rec3  = P_REC3.m();

	      HepLorentzVector P_REC4 = P_BEAM-Pi02.P4-Mu1.P4-Mu2.P4;
	      double m_rec4  = P_REC4.m();
	      
	      m_2pi0->column("rcm2mup1", m_rec3);
	      m_2pi0->column("rcm2mup2", m_rec4);
	      m_2pi0->column("m2pi2mu", mbc2);

	      double m2mupi1 = (Mu1.P4 + Mu2.P4 + Pi01.P4).m();
	      double m2mupi2 = (Mu1.P4 + Mu2.P4 + Pi02.P4).m();

	      m_2pi0->column("m2mupi1", m2mupi1);
	      m_2pi0->column("m2mupi2", m2mupi2);

	      m_2pi0->column("m2mu", P2DIMU.P4.m());

	      m_2pi0->dumpData();


	    
	  for(std::vector<HepParticle>::iterator it3 = myGamList.begin(); it3 != myGamList.end(); ++it3)
	    {

	      HepParticle& G = *it3;

              HepLorentzVector P_GREC = P_REC-G.P4;

	      HepLorentzVector P4_hbframe= G.P4;


	      double anglegwrthb = -10;
	      

	      //	      P4_hbframe.boost((P2PI.P4).boostVector());

	      P4_hbframe.boost(-P_REC.boostVector());



		  if (P4_hbframe.e()>=0.100&&P4_hbframe.e()<0.180){
		    if (my_ngam[0]<9){ my_ngam[0]=my_ngam[0]+1;}}
		  if (P4_hbframe.e()>=0.050&&P4_hbframe.e()<0.230){
		    if (my_ngam[1]<9){ my_ngam[1]=my_ngam[1]+1;}}
		  if (P4_hbframe.e()>=0.210&&P4_hbframe.e()<0.290){
		    if (my_ngam[2]<9){ my_ngam[2]=my_ngam[2]+1;}}
		  if (P4_hbframe.e()>=0.160&&P4_hbframe.e()<0.340){
		    if (my_ngam[3]<9){ my_ngam[3]=my_ngam[3]+1;}}
		  if (P4_hbframe.e()>=0.380&&P4_hbframe.e()<0.520){
		    if (my_ngam[4]<9){ my_ngam[4]=my_ngam[4]+1;}}
		  if (P4_hbframe.e()>=0.330&&P4_hbframe.e()<0.570){
		    if (my_ngam[5]<9){ my_ngam[5]=my_ngam[5]+1;}}
		  if (P4_hbframe.e()>=0.420&&P4_hbframe.e()<0.560){
		    if (my_ngam[6]<9){ my_ngam[6]=my_ngam[6]+1;}}
		  if (P4_hbframe.e()>=0.370&&P4_hbframe.e()<0.610){
		    if (my_ngam[7]<9){ my_ngam[7]=my_ngam[7]+1;}}
		  if (P4_hbframe.e()>=0.700&&P4_hbframe.e()<0.900){
		    if (my_ngam[8]<9){ my_ngam[8]=my_ngam[8]+1;}}
		  if (P4_hbframe.e()>=0.900){
		    if (my_ngam[9]<9){ my_ngam[9]=my_ngam[9]+1;}}



	      /*
	      std::cout << "P_BEAM=" << P_BEAM <<std::endl;
	      std::cout << "Pi1.P4=" << Pi1.P4 <<std::endl;
	      std::cout << "Pi2.P4=" << Pi2.P4 <<std::endl;
	      std::cout << "G.P4=" << G.P4 <<std::endl;
	      */


	      /*
      	      std::cout << "gammas=" << GamList.size() <<std::endl;
	      std::cout << "tracks=" << TrkList.size() <<std::endl;
	      std::cout << "pi0=" << pi0_count <<std::endl;
	      */


	      //	      if(G.P4.e() > 0.05 && GamList.size() < 20 && pi0_count < 10 && TrkList.size() < 20)
	      double m_grec = P_GREC.m();

	      if (m_grec < 9.2 || m_grec > 11.0) continue;

	      if (EvtMgr.count()>0) 
		{
		  m_2pi0g->column("expno",EvtMgr[0].ExpNo());
		  m_2pi0g->column("runno",EvtMgr[0].RunNo());
		  m_2pi0g->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
		}
	      m_2pi0g->column("ntr",     TrkList.size());
	      m_2pi0g->column("npi0",    Pi0List.size());

//	      m_2pi0g->column("ngam",    my_ngam_count);
	      m_2pi0g->column("ngam",    myGamList.size());
	      m_2pi0g->column("energy",E_BEAM);
	      m_2pi0g->column("evtr2",evt_R2);
              m_2pi0g->column("recm2pi", m_rec);
              m_2pi0g->column("m2pi",    P2PI.P4.m());

	      m_2pi0g->column("pxpi1",   Pi01.P4.px());
	      m_2pi0g->column("pypi1",   Pi01.P4.py());
	      m_2pi0g->column("pzpi1",   Pi01.P4.pz());
	      m_2pi0g->column("pxpi2",   Pi02.P4.px());
	      m_2pi0g->column("pypi2",   Pi02.P4.py());
	      m_2pi0g->column("pzpi2",   Pi02.P4.pz());


              m_2pi0g->column("recmpi1", m_rec1);
              m_2pi0g->column("recmpi2", m_rec2);

	      m_2pi0g->column("npostr",  Npos);
	      m_2pi0g->column("nnegtr",  Nneg);

	      m_2pi0g->column("mchrg", (P_CHG).m());

	      m_2pi0g->column("npipi",    my_npipi_count);

	      //	      m_new->column("eg",  G.P4.e());
	      //m_new->column("egboost",  P4_hbframe.e());

	      m_2pi0g->column("eg", G.P4.e());
	      m_2pi0g->column("egboost", P4_hbframe.e());
	      m_2pi0g->column("recmg2pi0", m_grec);

	      m_2pi0g->dumpData();



		}

	    



		}
	
	    }

	    
	  for(std::vector<HepParticle>::iterator it3 = myGamList.begin(); it3 != myGamList.end(); ++it3)
	    {

	      HepParticle& G = *it3;

              HepLorentzVector P_GREC = P_RECp1-G.P4;

	      HepLorentzVector P4_hbframe= G.P4;


	      double anglegwrthb = -10;
	      

	      //	      P4_hbframe.boost((P2PI.P4).boostVector());

	      P4_hbframe.boost(-P_RECp1.boostVector());



		  if (P4_hbframe.e()>=0.100&&P4_hbframe.e()<0.180){
		    if (my_ngam[0]<9){ my_ngam[0]=my_ngam[0]+1;}}
		  if (P4_hbframe.e()>=0.050&&P4_hbframe.e()<0.230){
		    if (my_ngam[1]<9){ my_ngam[1]=my_ngam[1]+1;}}
		  if (P4_hbframe.e()>=0.210&&P4_hbframe.e()<0.290){
		    if (my_ngam[2]<9){ my_ngam[2]=my_ngam[2]+1;}}
		  if (P4_hbframe.e()>=0.160&&P4_hbframe.e()<0.340){
		    if (my_ngam[3]<9){ my_ngam[3]=my_ngam[3]+1;}}
		  if (P4_hbframe.e()>=0.380&&P4_hbframe.e()<0.520){
		    if (my_ngam[4]<9){ my_ngam[4]=my_ngam[4]+1;}}
		  if (P4_hbframe.e()>=0.330&&P4_hbframe.e()<0.570){
		    if (my_ngam[5]<9){ my_ngam[5]=my_ngam[5]+1;}}
		  if (P4_hbframe.e()>=0.420&&P4_hbframe.e()<0.560){
		    if (my_ngam[6]<9){ my_ngam[6]=my_ngam[6]+1;}}
		  if (P4_hbframe.e()>=0.370&&P4_hbframe.e()<0.610){
		    if (my_ngam[7]<9){ my_ngam[7]=my_ngam[7]+1;}}
		  if (P4_hbframe.e()>=0.700&&P4_hbframe.e()<0.900){
		    if (my_ngam[8]<9){ my_ngam[8]=my_ngam[8]+1;}}
		  if (P4_hbframe.e()>=0.900){
		    if (my_ngam[9]<9){ my_ngam[9]=my_ngam[9]+1;}}



	      /*
	      std::cout << "P_BEAM=" << P_BEAM <<std::endl;
	      std::cout << "Pi1.P4=" << Pi1.P4 <<std::endl;
	      std::cout << "Pi2.P4=" << Pi2.P4 <<std::endl;
	      std::cout << "G.P4=" << G.P4 <<std::endl;
	      */


	      /*
      	      std::cout << "gammas=" << GamList.size() <<std::endl;
	      std::cout << "tracks=" << TrkList.size() <<std::endl;
	      std::cout << "pi0=" << pi0_count <<std::endl;
	      */


	      //	      if(G.P4.e() > 0.05 && GamList.size() < 20 && pi0_count < 10 && TrkList.size() < 20)
	      double m_grec = P_GREC.m();

	      if (m_grec < 9.2 || m_grec > 11.0) continue;

	      if (EvtMgr.count()>0) 
		{
		  m_pi0g->column("expno",EvtMgr[0].ExpNo());
		  m_pi0g->column("runno",EvtMgr[0].RunNo());
		  m_pi0g->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
		}
	      m_pi0g->column("ntr",     TrkList.size());
	      m_pi0g->column("npi0",    Pi0List.size());

//	      m_pi0g->column("ngam",    my_ngam_count);
	      m_pi0g->column("ngam",    myGamList.size());
	      m_pi0g->column("energy",E_BEAM);
	      m_pi0g->column("evtr2",evt_R2);
              m_pi0g->column("recmpi", m_rec1);

	      m_pi0g->column("pxpi1",   Pi01.P4.px());
	      m_pi0g->column("pypi1",   Pi01.P4.py());
	      m_pi0g->column("pzpi1",   Pi01.P4.pz());


	      m_pi0g->column("npostr",  Npos);
	      m_pi0g->column("nnegtr",  Nneg);

	      m_pi0g->column("mchrg", (P_CHG).m());


	      m_pi0g->column("eg", G.P4.e());
	      m_pi0g->column("egboost", P4_hbframe.e());
	      m_pi0g->column("recmgpi0", m_grec);

	      m_pi0g->dumpData();


	    }
    }
			}
		    }
		  tit4++;
		}
	      tit3++;

  eventcount++;
  if (eventcount % 100==0) printf("Event selection - %d/%d\n",savecount,eventcount); 
}

/*
 .......................................................................
 . FWJET2 - Subroutine to compute event sphericity 
 . 
 . Inputs    : NTRK -- number of tracks in momentum array 
 .           : PTRK -- array dimensioned (3,NTRK) containing 3-momenta 
 .           :         of tracks in the event 
 .           : EBEAM -- beam energy of the event 
 . Outputs   : H0 -- 0th moment of the event 
 .           : R1 -- H1/H0 of the event 
 .           : R2 -- H2/H0 of the event 
 .           : R3 -- H3/H0 of the event 
 .           : R4 -- H4/H0 of the event 
 . 
 . Calls     : LEGN(X) -- INTERNAL functions for Legendre polynomials 
 . Called    : Anyone who wants the Fox/Wolfram moments of the event 
 . Author    : F. Morrow  09/02/89  14.47.58 
 .             Transfer to C by kfjack 2/9,1999
 . 
 . Detailed description 
 .     FWJET2 determines the Fox/Wolfram moments of the event, up to 
 .     order 5. 
 .......................................................................
*/
void mumupi0pi0_hbook::fwjet2(int NTrk, float *PTrk, float EBeam, 
		     float *H0, float *R1, float *R2, float *R3, float *R4)
{
  int i, j;
  float sss ,cxang, ppari, pparj, delh0;
  float h, h1, h2, h3, h4;

  if (NTrk >= 2) {
    
    if (EBeam < 4.0 || EBeam > 6.5) EBeam = 5.29;
    sss = EBeam*EBeam*4;
    
    h=h1=h2=h3=h4=0;
    
    for (i=0;i<NTrk;i++) {
      ppari = sqrt(PTrk[i*3]*PTrk[i*3]+
		                PTrk[i*3+1]*PTrk[i*3+1]+
		   PTrk[i*3+2]*PTrk[i*3+2]);
      
      for (j=i;j<NTrk;j++) {
	
	if (i!=j) {
	  pparj = sqrt(PTrk[j*3]*PTrk[j*3]+
		                    PTrk[j*3+1]*PTrk[j*3+1]+
		       PTrk[j*3+2]*PTrk[j*3+2]);

	  cxang = (PTrk[i*3]*PTrk[j*3]+ 
		    PTrk[i*3+1]*PTrk[j*3+1]+
		   PTrk[i*3+2]*PTrk[j*3+2])/(ppari * pparj);
	  
	  delh0 = 2*ppari*pparj;
	  
	  h += delh0;
	  h1 += delh0 * cxang;
	  h2 += delh0 * LGNDR2(cxang);
	  h3 += delh0 * LGNDR3(cxang);
	  h4 += delh0 * LGNDR4(cxang);
	} else {
	  delh0 = ppari*ppari;
	  h += delh0;
	  h1 += delh0;
	  h2 += delh0;
	  h3 += delh0;
	  h4 += delh0;
	}
      }
    }
  }

  if (h > 1e-7f) {
    *H0 = h / sss;
    *R1 = h1 / h;
    *R2 = h2 / h;
    *R3 = h3 / h;
    *R4 = h4 / h;
  }else {
    *H0=*R1=*R2=*R3=*R4=-1.0;
  }
}


#if defined(BELLE_NAMESPACE)
} // namespace Belle
#endif

